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We consider an important class of self-assembly problems and using the formalism of stochastic 
thermodynamics, we derive a set of design principles for growing controlled assemblies far from 
equilibrium. The design principles constrain the set of structures that can be obtained under 
non-equilibrium conditions. Our central result provides intuition for how equilibrium self-assembly 
landscapes are modified under finite non-equilibrium drive. 


The fields of colloidal and nanoscale self-assembly have 
seen dramatic progress in the last few years. Indeed 
experimental and theoretical work has elucidated de¬ 
sign principles for the assembly of complex three dimen¬ 
sional structures m ■ Most of these advances, how¬ 
ever, are based on an equilibrium thermodynamic frame¬ 
work: the target structure minimizes a thermodynamic 
free energy [5]. Understanding the principles governing 
self-assembly and organization in far from equilibrium 
systems remains one of the central challenges of non¬ 
equilibrium statistical mechanics [SHE]. In this letter, 
we show that design principles can be derived for a broad 
class of non-equilibrium driven self-assembly processes. 
Our central result constrains the set of possible structures 
that can be achieved under a non-equilibrium drive. 

Imagine a self assembly process in which interactions 
amongst the various monomers are described by a set 
of energies E eq . The ratio of association and dissocia¬ 
tion rates is set by a combination of interaction ener¬ 
gies and chemical potentials {... /ij ... } of the monomers. 
This generic setup is sufficient to describe many self as¬ 
sembly processes. Examples include: growth of crys¬ 
tals from solution by nucleation [9], growth dynamics of 
cell walls pd], growth of multicomponent assemblies [4] 
and growth dynamics of biological polymers and fila¬ 
ments m■ The chemical potential controls the growth 
of the assembly. If the chemical potential is tuned to 
a coexistence value such that the assembly grows at an 
infinitesimally slow rate, then the structure of the assem¬ 
bly and its composition can be predicted by computing 
the equilibrium partition function and free energy G eq 
appropriate to the set of interaction energies. 

For values of the chemical potentials more favorable 
than the coexistence chemical potential, the assembly 
grows at a non-zero rate. In such instances, the struc¬ 
ture and composition of the growing assembly might not 
have sufficient time to relax to values characteristic of the 
equilibrium partition function 0EEMII. Defects are ac¬ 
cumulated as the self assembled structure grows at a non¬ 
zero rate. The time taken for a defect to anneal increases 
rapidly with distance from the interface of the growing 
structure. Due to the resulting kinetically-trapped states 
the crystal can assume structures very different from 
those representative of the equilibrium state [snsana. 

By applying the second law of thermodynamics and 
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FIG. 1. Schematic of the self-assembly problem considered 
in the paper. The probability of observing structural and 
compositional fluctuations in the growing assembly, p(u>) can 
be different from the canonical distribution p eq (u) specified 
by the interaction energies E eq . Our central result, Eq. |T] 
constrains the set of compositional fluctuations that can be 
achieved under a chemical potential drive 5p. 


the formalism of stochastic thermodynamics, we derive a 
surprising thermodynamic relation that is applicable to 
the above mentioned kinetic processes. This relation pro¬ 
vides constraints on the configurations that are achiev¬ 
able in a non-equilibrium self assembly process, 


d(N) t 

dt 


<5/i 


N 


> 0 , 


(1) 


where N is the size of the assembly at some instant 
of time, d(N) t /dt is the average rate of the growth of 
the assembly, pat(w) is probability distribution associated 
with a configuration ui in the growing assembly, p e ^(ui) = 
exp[— (E eq (u) — G^/ksT] is the equilibrium probabil¬ 
ity distribution obtained when the assembly is grown at 
an infinitesimally slow rate and D[p||(/] = f plnp/q > 0 
is the relative entropy between distributions p and q. 
The relative entropy is a measure of distinguishability 
between distributions p and q. It is zero only when the 
two distributions are identical and is nonzero otherwise. 
We have assumed that the chemical potential of each 
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monomer is the same and exceeds the equilibrium coexis¬ 
tence value by Sp (i.e. the concentrations of monomers in 
solution exceeds the equilibrium concentration required 
for assembly). The chemical potential difference Sp pro¬ 
vides the non-equilibrium driving force for self assembly. 

An alternative and presumably more practical formu¬ 
lation of the central result is in terms of interaction ener¬ 
gies. Imagine that we wish to generate compositions and 
structures in the growing assembly which are characteris¬ 
tic of a Hamiltonian E eS different from the Hamiltonian 
governing the interactions between species E eq . The cen¬ 
tral result places a bound on the minimum required ex¬ 
cess chemical potential Sp required to achieve such an 
assembly, 


d(N) t 

dt 


-G n + G% + (E eS - E eq ) N , f ' 

AT + 


> 0. 


( 2 ) 


Here Gat is the free energy of N particle system described 
by the Hamiltonian E eS . The free energies differences 
can either be computed directly from simulations or es¬ 
timated using an analytical framework. 

Hence, given a chemical potential drive Eq. [l] Eq. [2] 
constrain the set of allowed non-equilibrium structures 
found in the assembly. Alternately, given a target dis¬ 
tribution pn{w) or a target effective Hamiltonian E eS , 
the central result sets a bound on the minimum chemical 
potential driving force required to achieve the assembly. 
If the chemical potentials of the monomers can be var¬ 
ied, Sp in Eq. [l] and Eq. [2] is replaced by an average, 
(Ei i)n/N. The bound can be used to variationally 
optimize non-equilibrium driving forces {... pi ... } that 
result in configurations characteristic of a desired effec¬ 
tive energy landscape. 


THERMODYNAMIC BOUNDS FOR 
NON-EQUILIBRIUM SELF ASSEMBLY 


N and a microscopic configuration w at a particular in¬ 
stant of time. The entropy of the system S is given by 

s = ~k B Pt(u, N) In P t (w, N) (4) 

uj,N 


where ks is the Boltzmann constant. 

To proceed, we decompose the distribution P t (oj,N) 
as P t (uj,N) = P t (N)p N (oj) where Pt.{N) and Pn(u) are 
both normalized probability distributions. In perform¬ 
ing the decomposition we have assumed that the system 
has reached a steady state and that distribution of com¬ 
positional fluctuations for a given system size, and that 
Pn(uj) is independent of the time t. This decomposi¬ 
tion is the main assumption behind our theoretical ap¬ 
proach. With this assumption, we can associate an effec¬ 
tive energy function E eS (oj) with the distribution p^( w), 
E eS (uj) ~ — lnpAr(o;). This is simply a statistical defin¬ 
ing relation for the energy function (analogous to a po¬ 
tential of mean force). The energy function E eS (ui) does 
not control the dynamics of the system. This effective en¬ 
ergy function can have a form very different from that of 
the interactions specified by the interaction Hamiltonian 
E eq (uS). 

Using this decomposition and the total average entropy 
production rate in a self assembly process, dStot/dt can 
be written as 
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A closer investigation reveals that (eaiss) > 0[? ]. 

Hence, in order to have a growing assembly with a dis¬ 
tribution of compositional fluctuations characteristic of 
the energy function E eS , the chemical potential differ¬ 
ence Sp must at least exceed (edi ss )> 


We now outline the derivation of the central result. 
This derivation is based on a generalization of the frame¬ 
work introduced in Ref. m ■ We begin by invoking 
one of the central results of stochastic thermodynam¬ 
ics Haora for a system in contact with a heat bath, 

d utnt dS dS P / \ 

< 3) 

where S is the entropy of the system and S e is the en¬ 
tropy due to energy exchange between the system and 
environment. Eq. [3] is simply a statement of the second 
law of thermodynamics. 

The growing self assembling structure will be the ther¬ 
modynamic system for the purposes of this paper. The 
system is allowed to exchange material and energy with 
the environment. Let Pt(u} 7 N) denote the probability 
distribution associated with the observing a system size 


dp A (^diss) • (6) 

A simple algebraic manipulation allows us to rewrite 
Eq. [b] as Eq. [l] When the chemical potentials of the 
monomers contributing to the assembly are unequal, Sp 
is replaced by the average excess chemical potential for 
the growing assembly, (JA Spi )n/N. 

The structure of Eq. [6] is reminiscent of the varia¬ 
tional Gibbs-Bogoliubov-Feynman inequality generalized 
for non-equilibrium systems. As such it can be applied to 
obtain bounds on optimal non-equilibrium driving forces 
that can maximize the yield of desired structures. Fi¬ 
nally, our central results show how the average entropy 
production rate in a self assembly process is coupled to 
the compositional fluctuations in the self assembled sys¬ 
tem. Recent work by Barato, Seifert and coworkers has 
demonstrated that the entropy production rate bounds 
the fluctuations of various dynamical currents generated 
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FIG. 2. Application of the bounds to ID fiber growth model (a) Schematic of the ID polymer growth process. The polymer is 
assembled from a bath composed of two monomer types (Red and Blue), (b) Schematic of the effective Markov state model. 
This Markov state model resolves the nature of terminal bond in the self-assembled system (vertical rungs), S = bond between 
like monomers and D = bond between unlike monomers, and the number of particles in the self assembled system (horizontal 
axis). The rates of transitions between the mesoscopic states resolved in this effective model depend on the composition of the 
assembly and are described in the (SI). We demonstrate that this effective model captures the dynamics and compositional 
fluctuations of the self assembly process, (c) Comparison between the lower bounds of 5[i obtained from Eq. [fijand Eq. [ 7 ] and 
the value of 5[i obtained from simulations. The equilibrium domain length, £ 0 , for this plot is 50. The red shaded region 
represents the region that is disallowed by Eq. [7] Our thermodynamic bounds are valid and inspite of their minimal nature, 
do reasonable job of predicting the compositional fluctuations in the assembly. 


in a non-equilibrium processes [20M22] . In the next sec¬ 
tion, using a model system, we demonstrate that the 
bound imposed by the entropy production rate on cunru- 
lants of fluctuations of the growth rate of the assembly 
lead to a set of tighter design principles, 

5n - (e diss ) > , (7) 

- t ^ §N ^ 


where (N) = N/t denotes the average rate of growth of 
the assembly, (( SN ) 2 ) denotes the variance in the growth 
rate of the assembly. 

We will now illustrate the result on two space filling lat¬ 
tice based model systems. In each example, the distribu¬ 
tion of compositional fluctuations in the growing assem¬ 
bly is modified when the assembly is grown at a non-zero 
rate. Our bound acts as a design principle and provides 
intuition into how the equilibrium landscape of compo¬ 
sitional fluctuations is modified under a non-equilibrium 
drive. We emphasize that while these examples focus on 
distributions of compositional fluctuations, our central 
results are valid more generally and also apply to statis¬ 
tics of structural fluctuations [9]. We choose the lattice 
based models as examples because they are analytically 
tractable and serve to illustrate the utility of Eq. [6] and 
Eq.0 


EXACTLY SOLVED MODEL OF 
NON-EQUILIBRIUM POLYMER ASSEMBLY 

As our first example, we will consider a model of 
one-dimensional polymer assembly. For simplicity, the 
system that we are studying will contains two types 
of monomers with equal concentrations. In addition, 
only nearest neighbors interaction is allowed between 
monomers. In scenarios where the polymer is nucleated 
from a solution, the concentration of monomers or their 
chemical potential provides the non-equilibrium driving 
force for the assembly. If the average length of monomer 
domains in a polymer at equilibrium is Co, our results 
impose a constraint on the excess chemical potential re¬ 
quired to achieve assemblies in which the average length 
of monomer domains is £ (SI), 


5/j, > 


(1 



£-1 
Co - 1 



— (cdiss) 


( 8 ) 


where 5/i is the excess chemical potential and we have 
assumed that the two monomers have the same chem¬ 
ical potential. The compositions in such an assembly 
differ from those obtained close to equilibrium ~ 0) 
when the polymer grows at an infinitesimally slow rate. 
Our thermodynamic bound doesn’t depend on the ki¬ 
netic model chosen. Further, as discussed previously, in 
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cases with unequal chemical potentials, Eq.[8]can be gen¬ 
eralized by simply replacing Sfj with the average excess 
chemical potential. 

We theoretically verify the bound for the polymer as¬ 
sembly problem by mapping the dynamics of the self as¬ 
sembly process onto an effective Markov state model (Fig 
[2] (b)). This effective model resolves the terminal bond in 
the self-assembled system (vertical rungs) and the num¬ 
ber of particles in the self assembled system (horizontal 
axis). The details of the various transition rates between 
the mesoscopic states are provide in the supplemental in¬ 
formation (SI). In the SI, we also demonstrate that the 
Markov state caricature accurately describes the dynam¬ 
ics and compositional fluctuations of the self assembly 
process. Further, we demonstrate that the expression 
for the entropy production rate, dS/ot /dt computed from 
this Markov state caricature indeed equals 


d*Stot 

dt 


dN 

dt 


(5/J, (fidiss) ) 


(9) 


where (edi ss ) is defined in Eq. [8] This calculation verifies 
the thermodynamic bound and its underlying assump¬ 
tions for the one dimensional polymer assembly model. 

The Markov state caricature allows us to map the self 
assembly process onto the dynamics of a random walker 
on a graph. In the SI, we verify that the effective mean 
field Markov state caricature manages to reasonably de¬ 
scribe even the second cumulant of fluctuations in the 
growth rate. Barato, Seifert and coworkers have recently 
demonstrated pOH22] that fluctuations of various dynam¬ 
ical currents in such Markov state process are bound by 
the entropy production rate. As demonstrated in the SI, 
we adapt these bounds to show that 


Sfi > 


(1 



£o-l 



2(N) 

t((SN) 2 ) 


( 10 ) 


where (N) is the average growth rate of the polymer and 
(( SN ) 2 ) is the variance of the growth rate. 

In Fig. ^ c), we demonstrate the two thermodynamic 
bounds in Eq. [8j Eq. [lO] for a particular kinetic model of 
polymer growth 12.11 . Including information from higher 
order cumulants does indeed strengthen the bound. 
These demonstrations prove the theoretical validity and 
effectiveness of our bound and pave the way for applica¬ 
tions to more complex systems. 

Specifically, for systems with more than two monomer 
types or for monomers with longer ranged interactions, 
it is more convenient to express (ediss) i n terms of the 
actual E eq and effective E e g interaction Hamiltonians. 
While effective Markov state caricatures for such sys¬ 
tems can be complex and analytically intractable, our 
thermodynamic bounds are easily accessible and provide 
very useful intuition for the compositional fluctuations 
of the non-equilibrium assembly. We now demonstrate 
these ideas by considering a two-dimensional analogue of 
the fiber growth process described above [M] (Fig. [if] (a)). 


NON-EQUILIBRIUM TWO DIMENSIONAL SELF 
ASSEMBLY PROCESS 

As a natural generalization of the one-dimensional 
polymer model, we imagine a two-dimensional assem¬ 
bly growing at one end as described in Fig [3] (a). We 
again imagine two monomer types with interactions be¬ 
tween like monomers set by the energy scale e s and 
interactions between unlike monomers set by the en¬ 
ergy scale ed- The details of the kinetic model are de¬ 
scribed in Methods. When the chemical potential is set 
to the intensive equilibrium free energy of the assembly, 
gcoex ~ 2e s — ksT In 2, the assembly doesn’t grow on av¬ 
erage. Assigning spin s = 1 to the red monomers and 
spin s = — 1 to the blue monomers, the statistics of the 
compositional fluctuations in the equilibrium assembly 
are equivalent to that of an Ising model with coupling 
constant J = e % ed . We will work in regimes where the 
coupling constant is above the critical coupling constant 
of the Ising model. J > J c ~ 1/2.26. We are interested in 
the statistics of compositional fluctuations as the growth 
rate of the assembly is tuned by changing the chemical 
potential. 

We begin by assuming that the compositional fluctu¬ 
ations in the growing assembly can be described by an 
effective nearest neighbor coupling constant J e g ■ We per¬ 
formed simulations in which we extracted the value of 
Jeff as a function of the excess chemical potential Sfj by 
analyzing compositional fluctuations in the growing as¬ 
sembly (Methods). The values of J e g obtained from sim¬ 
ulations are plotted in Fig. [3])b). The underlying Ising 
model had a coupling constant J = 0.75. When the ex¬ 
cess chemical potential is close to Sfj ~ 0.30, the value 
of the effective coupling constant computed from simu¬ 
lations reaches the critical value of the two dimensional 
Ising model, J c . 

Using Eq.[6]and Eq.[7]we analytically computed a lower 
bound Jb on the values of the effective coupling constant 
as a function of Sfj,. In Fig. [ 3 ] (b), we plot the difference 
J e ff — Jb for various values of the excess chemical poten¬ 
tial. The bound obtained from Eq. [7] does a reasonable 
job of predicting the effective coupling constant. Includ¬ 
ing information about higher order cumulants again im¬ 
proves the bound and provides a close to accurate esti¬ 
mate of the effective coupling constant J e g. In partic¬ 
ular, the bound obtained by including information form 
the higher cumulants (Eq. [7]) does particularly well for 
values of Sfj around the critical value, 5fi c ~ 0.30. This 
close agreement between the bound predicted by theory 
and the effective value computed from simulations in the 
critical regime is surprising given the minimal nature of 
our theory. 

We found similar results for other values of the cou¬ 
pling constant J. Strikingly, for all the regimes consid¬ 
ered, the bounds are particularly accurate when the as¬ 
sembly is close to criticality. The results illustrate the 
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FIG. 3. Application of the bounds to the non-equilibrium two dimensional growth process (a) We assemble a two dimensional 
structure from a monomer bath containing red and blue blocks. The energy of interaction between similar monomers is e s while 
that between dissimilar monomers is ej. As discussed in the text, when the assembly is grown at equilibrium, the statistics of 
compositional fluctuations in the assembly are equivalent to that of an Ising magnet with coupling constant J = . (b) 

For assemblies grown away from equilibrium, we analyzed the compositional fluctuations and using the analytical solution for 
a two dimensional Ising magnet, compute the effective coupling constant J e g as a function of the excess chemical potential Sfj,. 
The equilibrium coupling constant was set to J = 0.75. As the system deviates more from equilibrium by increasing Sfj, the 
J e ff decreases and reaches a critical value J c for Sfj > 0.30 . (c) We computed bounds for the effective coupling constant, Jj, 
using Eq. |6]and Eq. [ 7 ] at J = 0.75. As evidenced by the difference J e s — Jb, the bounds are able to offer reasonable intuition 
for the statistics of compositional fluctuations in the growing non-equilibrium assembly. The bounds work particularly well for 
values of excess chemical potential Sfi close to the critical value 5fi c « 0.30. This result is striking given the minimal nature of 
our model. 


utility of our central relations as design principles for 
non-equilibrium self assembly. 

APPLICATIONS AND CONCLUSION 

Using the framework of stochastic thermodynamics 
and physically justified approximations, we have put for¬ 
ward a general predictive framework for non-equilibrium 
self assembly. Our central results, Eq. [6] and Eq. [7] are 
built around a simple and intuitive expression for en¬ 
tropy production in non-equilibrium self assembly pro¬ 
cess. This expression relates compositional fluctuations 
in the self assembled system to the chemical potential 
excess driving assembly process. Given a set of pro¬ 
grammed interactions between monomers and structures 
specified by the energy function E eq , our central result 
constrains the set of allowed structures that can be ob¬ 
tained under a non-equilibrium drive. A unique feature 
of our result is its thermodynamic character: the bound 
doesn’t depend on the details of the kinetic model cho¬ 
sen. Many recent theoretical and experimental studies 
have postulated that non-equilibrium protocols can po¬ 
tentially be used to overcome bottlenecks and enhance 
the yield of desired self assembled structures [25] . Our 
framework is ideally suited to verify these postulates and 
explore connections between energy consumption and im¬ 
proved yield of target structures. 

We anticipate that our results will find important ap¬ 
plications in studies of crystal nucleation m, self as¬ 
sembly of metal organic frameworks 126 1, non-equilibrium 


roughening transitions, and the synthesis of nano crystals 
using DNA and other biopolymers as building materials. 
In all these instances, the assembly process can occur 
far from equilibrium depending on the effective chemical 
potentials imposed. Our result provides a bound on the 
chemical potentials under which the desired structures 
can be robustly obtained even far from equilibrium. It 
also provides crucial intuition into how the equilibrium 
landscape Eng and the capacity of the system to store 
structures m is modified under non-equilibrium condi¬ 
tions. 

Finally, the variational structure of our central results 
provides a framework to optimally tune the chemical po¬ 
tential driving forces so as to maximize the yield of de¬ 
sired structures. In this context, our results complement 
recently derived theoretical results that suggest a het- 
eregeous set of excess chemical potentials is crucial for 
robust nucleation of large complex nanoscale assemblies 
using DNA origami [2B„ 22] • 

A mismatch of chemical potentials provides the non¬ 
equilibrium driving force in Eq. Our results also 
apply to non-equilibrium error correction mechanisms 
for biopolymer assembly in which the assembly pro¬ 
cess is driven by consumption of energy rich ATP/GTP 
molecules. In such instances, the Sfj term is replaced 
by the appropriate average energy consumption rate. In 
these biophysical applications, the non-equilibrium driv¬ 
ing force generates structures that posses correlations 
greater than those in equilibrium. While such fibers are 
not achievable with the kinetic protocols described above, 
the thermodynamic nature of our bounds allows their 
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application to these biophysical error correction mecha¬ 
nisms. Indeed, Sartori et al [30l have recently derived a 
similar bound for a specific model of kinetic proofread¬ 
ing. Our general thermodynamic bounds, particular that 
in Eq. [ 7 ] can elucidate the tradeoffs between speed ac¬ 
curacy and energy consumption in such non-equilibrium 
biochemical replication processes. 

METHODS 

In our ID polymer growth simulation, we consider 
a polymer growing in a bath containing two types of 
monomer. The concentrations of both types are equal 
to c. The chemical potential f.i is related to the concen¬ 
tration c by the relation: /it = —/csTln(c). The nearest 
neighbor monomers are allowed to interact with energy 
e s if they are the same and with energy if they are 
different. We perform kinetic Monte-Carlo simulations 
in which monomers are added to the polymer with the 
rate c (or exp(—/?//)). Monomers are allowed to dissoci¬ 
ate from the polymer with the rate exp(—/Je,;) (i here can 
be either s or d) which depends solely on the composition 
of the fiber. /3 is set to 1. We use </> to denote the proba¬ 
bility that a particular bond in the assembly is between 
two like monomers. It is related to the correlation length 

£ by: 0=1- 1/C- 

The simulations for the non-equilibrium 2D assembly 
problem follow the same rules as in the ID case. However, 
we used the standard Metropolis Monte-Carlo algorithms 
for the 2D case. For the purposes of computing the ef¬ 
fective coupling constant J e g we labeled the red particles 
as spin up s = 1 and blue particles as spin down s = — 1. 
We then computed the average per unit magnetization m 
in the growing assembly and used the Onsager solution 
to extract the value of J e g. This procedure is only valid 
for J gff V J c . 
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